%% EXERCISE 33: Infinite plate with a hole under uniaxial stress
clear all
close all

format short
addpath('ex33_toolbox_local')

disp('================ EXERCISE 33: Infinite plate ===================')

%% Setting parameters and options
a     = 1;
L     = 5;
E     = 1e1;            %plane stress
nu    = 0.25;         %plane stress
gamma = 1.6;

E      = E / (1.0-nu*nu);        %plane strain
nu     = nu / (1.0-nu);          %plane strain

lambda = nu*E / (1.0-nu*nu);
kappa  = (3.0-nu)/(1.0+nu);
mu     = 0.5 * E / (1.0+nu);


parameters.a      = a;
parameters.L      = L;
parameters.E      = E;
parameters.nu     = nu;
parameters.gamma  = gamma;
parameters.lambda = lambda;
parameters.kappa  = kappa;
parameters.mu     = mu;

Naxis = [ 6 12 24 36 48 72 96 192];
N = Naxis(2);

% Nodal Points
x_nodes = ex33_mesh_generation2(N,a,L);
nPts    = length(x_nodes);


%% Analytical computations

% stress and displacement
s_xx_ana = zeros(nPts,1);
s_yy_ana = zeros(nPts,1);
s_xy_ana = zeros(nPts,1);
u_ana    = zeros(nPts,2);
for i=1:nPts
  x = x_nodes(i,1);
  y = x_nodes(i,2);
  [s_xx_i,s_yy_i,s_xy_i] = ex33_AnalyticalStress(x,y,a);
  s_xx_ana(i) = s_xx_i;
  s_yy_ana(i) = s_yy_i;
  s_xy_ana(i) = s_xy_i;
  [u_x u_y]   = ex33_AnalyticalDisplacement(x,y,a,mu,kappa);
  u_ana(i,:)  = [u_x u_y];
end

figure(1)
plot(x_nodes(:,1),x_nodes(:,2),'sb')
axis equal
hold on
plot(u_ana(:,1)+x_nodes(:,1),u_ana(:,2)+x_nodes(:,2),'.r')
hold on
quiver(x_nodes(:,1),x_nodes(:,2),u_ana(:,1),u_ana(:,2),'k')
legend('Reference','Deformed','Displacement field')
title('Analytical Deformation mapping')

% traction
t_ana          = [];
parameters.bd1 = [];
parameters.bd2 = [];
parameters.bd3 = [];
parameters.bd4 = [];
for i=1:nPts
  [t_ t_bd parameters] = ex33_AnalyticalTraction(x_nodes(i,1),x_nodes(i,2),s_xx_ana(i),s_yy_ana(i),s_xy_ana(i),parameters,i);
  if t_bd == 1
    t_ana = [t_ana; t_];
  end
end
id_bd = union(parameters.bd1,union(parameters.bd2,union(parameters.bd3,parameters.bd4)));


figure(2)
quiver(x_nodes(id_bd,1),x_nodes(id_bd,2),t_ana(:,1),t_ana(:,2),'r')
title('Traction field')
hold on
axis equal


%% Generation of Gauss points for numerical integration
%cubature order: 1 (1 GPts), 2 (3 GPts), 4 (6 GPts), 5 (7 GPts), 6 (12 GPts)
optSample.pruning = 0;
optSample.orderGL = 2;
optSample.quality = 0.1;
[x_samples w_samples] = MakeGLSamples2D(x_nodes, optSample);
eps   = 1.0e-7;

sPts = size(x_samples,1);
k=1;
for i=1:sPts
  x = x_samples(i,1);
  y = x_samples(i,2);
  r = sqrt(x*x+y*y);
  if r<a %to avoid add material in the hole
    %disp('x_samples is not considered');
  else
    x_samples(k,:) = x_samples(i,:);
    w_samples(k,:) = w_samples(i,:);
    k=k+1;
  end
end
sPts = k-1;
x_samples(k:end,:)=[];
w_samples(k:end)=[];

figure(3)
plot(x_nodes(:,1),x_nodes(:,2),'ob')
hold on
plot(x_samples(:,1),x_samples(:,2),'xk')
legend('Nodes','Samples')
title('Nodal and sample points')
axis equal


%% Computation of shape functions

optLME.dim   = 2;
optLME.verb  = 0; %0:off
optLME.TolNR = 1.e-14;
optLME.knn   = 1;
optLME.Tol0  = 1.e-10;
optLME.gamma = gamma;      % value of gamma to compute LME
[beta_n range_n h_n] = NodalThermalization(x_nodes, optLME, id_bd);

% h       = 0.25;
% beta_n  = ones(length(x_nodes),1)'*optLME.gamma/h/h;
% range_n = sqrt(-log(optLME.Tol0)./beta_n)';
optLME.beta_n = beta_n;
optLME.range_n= range_n;
[p_lme, dp_lme, s_nears]= ex33_wrapper_lme(x_nodes,x_samples,optLME);

clear hp_lme;

u_num = ex33_SolveSystemHolePlate(dp_lme,s_nears,x_nodes,x_samples,w_samples,parameters,optLME);
u_num = reshape(u_num,2,length(x_nodes))';

figure(4)
plot(x_nodes(:,1),x_nodes(:,2),'.k')
hold on
plot(u_ana(:,1)+x_nodes(:,1),u_ana(:,2)+x_nodes(:,2),'.r')
hold on
plot(u_num(:,1)+x_nodes(:,1),u_num(:,2)+x_nodes(:,2),'.b')
hold on
title('Deformation Mapping (nodal points)')
legend('Reference','Analytical','Numerical');
axis equal


figure(5)
plot(u_ana(:,1),u_ana(:,2),'.r')
hold on
plot(u_num(:,1),u_num(:,2),'.b')
hold on
title('Displacement (nodal points)')
legend('Analytical','Numerical');
axis equal

u_num_s = SamplesSolution(u_num,s_nears,p_lme);

u_ana_s = zeros(sPts,2);
for i=1:sPts
  x = x_samples(i,1);
  y = x_samples(i,2);
  [u_x u_y]    = ex33_AnalyticalDisplacement(x,y,parameters.a,parameters.mu,parameters.kappa);
  u_ana_s(i,:) = [ u_x u_y];
end
L2 = sqrt(sum(sum((u_num_s-u_ana_s).^2,2).*w_samples))/sqrt(sum(sum(u_ana_s.^2,2).*w_samples));
fprintf(1,'L2-norm error = %8.2e\n',L2);

figure(6)
plot(u_ana_s(:,1)+x_samples(:,1),u_ana_s(:,2)+x_samples(:,2),'.r')
hold on
plot(u_num_s(:,1)+x_samples(:,1),u_num_s(:,2)+x_samples(:,2),'.b')
hold on
title('Deformation Mapping (sample points)')
legend('Analytical','Numerical');
axis equal

